Self-interaction correction in multiple scattering theory 
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We propose a simplified version of self-interaction corrected local spin-density (SIC-LSD) approx- 
imation, based on multiple scattering theory, which implements self-interaction correction locally, 
within the KKR method. The multiple scattering aspect of this new SIC-LSD method allows for the 
description of crystal potentials which vary from site to site in a random fashion and the calculation 
of physical quantities averaged over ensembles of such potentials using the coherent potential ap- 
proximation (CPA). This facilitates applications of the SIC to alloys and pseudoalloys which could 
describe disordered local moment systems, as well as intermediate valences. As a demonstration of 
the method, we study the well-known a-7 phase transition in Ce, where we also explain how SIC 
operates in terms of multiple scattering theory. 
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I. INTRODUCTION 

The self-interaction corrected local spin density (SIC- 
LSD) approximation^'^ has proved to be a useful scheme 
to describe static correlations in strongly correlated elec- 
tron systems. In particular, it can determine whether an 
electron is delocalized or localized, i.e., whether its or- 
bital is part of the valence states or not. This leads to 
a determination of the number of valence states and a 
nominal valence, as demonstrated by numerous calcula- 
tions on rare earths, actinides, transition metal oxides, 
including the parent compounds of the high Tc materials 
and the CMR materials.^''^-^-^-'^-®-^ 

The full SIC-LSD scheme is unfortunately difficult to 
implement.^P This is due to the repeated transformations 
from reciprocal space (k-space) to real space to evaluate 
the self-interaction potential and the back transforma- 
tions to k-space to solve the band structure problem. So 
far most applications of the full SIC-LSD formalism have 
been implemented in the LMTO-ASA (linearized muffin- 
tin orbitals) band structure methodic In this paper a 
simpler but more versatile scheme is developed and im- 
plemented within multiple scattering theory, in the Kor- 
ringa, Kohn, Rostocker (KKR) formulation. Its main ad- 
vantage, thanks to a straightforward determination of the 
Green's function, is a possible generalization to alloys via 
the coherent potential approximation (CPA) J^ii2iiiiiSii& 
Since a single-site approximation underpins this new for- 
mulation, in what follows it is referred to as a local self- 
interaction correction (LSIC) formalism. It is based on 
the experience with the full SIC-LSD implementation 
showing that to better than 98% the electron is local- 
ized on the site under consideration, which justifies the 
single-site approximation. Whilst in the full LMTO-ASA 
implementation, the representation of the localized or- 



bitals over a real-space cluster determines the extent of 
these orbitals, in the present scheme the degree of local- 
ization is determined by the energy dependence of the 
single-site phase shift, in particular the width of its reso- 
nance, corresponding to the localized electron. A broader 
resonance would imply reduced localization. 

The SIC-LSD formalism for solids has been developed 
into a scheme that treats both localized and delocalized 
electrons on equal footing. The decision whether an elec- 
tron should be considered as localized or delocalized is 
based on a delicate balance between the energy gain due 
to the inclusion of the self-interaction correction (local- 
ization) energy and the energy loss in band or hybridiza- 
tion energy.'^ While this methodology has been success- 
ful in differentiating localized from delocalized electrons, 
i.e., a dual character of the electron, it does not describe 
the interesting crossover between localized and delocal- 
ized states which occurs for example in heavy fermion 
systems. Our aim here is to develop a theory which 
describes local fluctuations of the electronic configura- 
tions between that where an electron can be said to be 
localized and another where an electron is delocalized. 
It will be shown that the present local formulation of 
SIC-LSD readily lends itself to be the basic idea of such 
a development. The origin of our approach goes back 
to the invention and use of the coherent potential ap- 
proximation to describe the charge and spin fluctuations 
about the Hartree-Fock solution of the Hubbard model 
by Hubbard himself (Hubbard III approximation) .4Li2iiS 
The present implementation of this idea rests on its gen- 
eralization to account for the corresponding fluctuations 
about the local density approximation (LDA) to the first- 
principles density functional theory (DFT)i ^°i^^i^^ In the 
literature this generalization, when applied to spin fluc- 
tuations, is referred to as the KKR-CPA implementa- 



tion of the disordered local moment (DLM) pictureiS2i24 
Indeed, the present work can be considered as the fur- 
ther elaboration of this basic idea in which the LSIC 
replaces the LDA as the local description of the elec- 
tronic structure and the attention is being focused on 
the valance fluctuations. Interestingly, it is now well es- 
tablished that Hubbard's so-called "alloy analogy" ap- 
proximation, which prompted the use of the CPA, leaves 
out of consideration some very important fluctuations. 
The most significant of these are those which give rise to 
a Kondo-like resonance at the Fermi energy in the case 
of the Hubbard model and correspond to such a qualita- 
tively new physics as the Mott transition. The relevance 
of this in the present context is that such fluctuations 
are well described by the dynamical mean field theory 
(DMFT)^'"^ whose static limit, for the Hubbard model, 
is precisely the "alloy analogy" CPA of Hubbard. -^^ This 
point was particularly clearly explained in the recent pa- 
per of Kakahashi.^® Consequently, it is reasonable to re- 
gard our LSIC based KKR-CPA-DLM calculations as in- 
vestigations of the static limit of a yet undeveloped first 
principles DMFT. In what follows when we refer to the 
need to include dynamical effects in the theory it is the 
above theoretical considerations we will have in mind. 

The paper is organized as follows. In Sec. [H] we out- 
line the physical picture underlying the present approach. 
In Sec. mil a general formulation of SIC-LSD, following 
Perdew and Zunger,-'^ is briefly summarized with refer- 
ence to some aspects of the LMTO-ASA implementation 
based on the Wannier function representation of local- 
ized orbitals.^*' In Sec. IIVI the formalism of the local 
self-interaction corrected local spin density (LSIC-LSD) 
within multiple scattering theory is described in detail. 
There we concentrate on the phase shifts and single-site 
Green's function from which the SIC charge and poten- 
tial, corresponding to localized electron states, are cal- 
culated within the KKR method. Since the latter can 
be easily extended to include coherent potential approx- 
imation, section briefly summarizes its most impor- 
tant equations in terms of the multiple scattering quan- 
tities. In section IVII the formalism is extended to finite 
temperatures. The potential and versatility of the LSIC 
method is demonstrated on the application to the a — 7 
phase transition in Ce. In section IVllI we first discuss 
the f-phase shifts, total energies, lattice parameters, den- 
sities of states (DOS), and spectral functions at T=OK 
for the a— and 7— phases, as obtained from the LSIC- 
KKR method. Wherever appropriate we compare with 
the results of the full SIC-LSD implementations within 
LMTO-ASA.^" In this section we also present calcula- 
tions for finite temperatures and the full phase diagram 
of the a — 7 phase transition. Both the CPA and DLM 
are utilized to accomplish the latter, and to illustrate 
how the present approach is capable to describe both spin 
and valence fluctuations at finite temperatures. Section 
IVIIIl is devoted to various aspects of the present approach 
and among them a consideration of how intermediate va- 
lence could be realized within the present implementa- 



tion, which motivates a possible generalization to include 
dynamics, as outlined in section HXI The paper is sum- 
marized in section Kl 



II. PHYSICAL PICTURE 

In the present formulation of SIC, we adopt the phys- 
ical picture of multiple-scattering theory, where a solid 
is represented by an array of non-overlapping scatter- 
ing centers. The electronic motion is then broken down 
into a sequence of scattering events and a free propaga- 
tion in between. The most useful concept of this method 
is a phase shift, describing scattering of electrons from 
ions, the scattering centers in a solid. If a phase shift 
is resonant it is reminiscent of a bound state at positive 
energies, i.e., above the zero of the potential which in 
our case is the mufSn-tin zero. The energy derivative of 
the phase shift is related to the Wigner delay-time. If 
this is large the electron will spend a long time on the 
site. Such 'slow' electrons will be much more affected by 
the spurious self-interaction and therefore should see a 
Sl-corrected potential. 

In most systems, where the electrons are truly delocal- 
ized, the self-interaction contribution to the potential is 
negligible and therefore the LDA is an excellent approx- 
imation. When the phase shift has a resonance one has 
to calculate the self-interaction correction for this (/,to) 
angular momentum channel. This is accomplished by 
calculating the one-electron charge density for this chan- 
nel, defining the charge density for the self-interaction 
correction. From this one can readily calculate the self- 
interaction potential which has to be added to the LSD 
potential, and then the new phase shifts are calculated 
for the total (SIC-LSD) potential. This has to be carried 
out m channel by m channel for a given angular momen- 
tum I. The minimum of the total energy will determine 
the optimum configuration of {l,m) channels to be self- 
interaction corrected. Therefore to each of the m chan- 
nels one can assign two potential functions, V^g (r) 
and VgJ'^(r). The formalism determining the energy 



functional associated with the potential V^jj 



SIC-LSD 



(r) is 



briefly outlined in the next section. It should be men- 
tioned here that if the total energies corresponding to 
these two different potentials are sufficiently close, one 
can envisage dynamical effects to play an important role 
as a consequence of possible tunneling between these 
states. We shall return to this point in the later sections. 



III. SIC-LSD FORMALISM 

It has been pointed out by Perdew and Zunger^ that 
density functional theory (DFT) schemes, like the local 
spin density approximation, suffer from a spurious self- 
interaction of the electrons with themselves. In princi- 
ple, this self-interaction term should vanish exactly, as it 
does in the Hartree-Fock theory. In practice, however. 



this cancellation is incomplete. Perdew and Zunger sug- 
gested an approximate solution to this problem, which 
was constructed for finite systems but is here extended 
to solids in a different way as compared to previous im- 
plementations for solids J 

The usual representation of the total energy within the 
LSD-DFT formalism in the Kohn-Sham approach^^ is 

occ 

Oca 

+Eu[n]+E^!^[nT,n^], (1) 

where ^ao-'s are the Kohn-Sham orbitals, (acr is a com- 
bined index labeling the orbital and spin ("f or |), respec- 
tively), Uaa = l^acrP, "-cr = 1]°'''' "-qct, U ^ 11^ + n| . E'cxt 

is the external field energy functional, £'h is the Hartree 
energy 



Eii[n] = fd^r fd^r 



3. Ls / n{r)n{r') 



(2) 



and E^^^ is the LSD approximation to the exchange- 
correlation energy functional. On the basis of the above, 
Perdew and Zunger proposed a self-interaction corrected 
LSD on an orbital by orbital basis 

occ 
- ^(^HKa] + Sj^i^^Ka, 0]),(3) 

by subtracting explicitly the self-Coulomb and self- 
exchange and self-correlation energy of all occupied or- 
bitals. This correction restores the property that the 
true functional E[n] should have, namely that the self- 
Coulomb energy exactly cancels the self-exchange and 
self-correlation energy for every single orbital. En [n^cr] + 
-^xc'''^* ["acr , 0] = 0. This leads to an orbital dependent 
SIC-potential seen by an electron in orbital (j>aa, 



ySIC-LSD 
eff,aCT 



(r) = K=,t(r)-HyH[n](r)-fl/c. K,nx](r) 

^^ v ' 

- T4jK.](r) - K'^e'^K., 0](r), (4) 
^^ V ' 

with the external lattice potential Kixtlr), and 

, n(r') 



FHN(r) = 2 /"dV 



14^-K,nx](r) 



|r — r' 



(5) 
(6) 



This self-interaction correction vanishes exactly only for 
extended states. In order to apply the SIC scheme to 
solids, the approach by Perdew and Zunger has to be 
generalized. This involves simultaneously a Wannier rep- 
resentation of the orbitals, necessary to determine Uaa of 



Eq. Q, and a Bloch representation to solve the band 
structure problem. Furthermore, the Wannier functions 
are required to fulfill the localization criterion which en- 
sures that the energy functional is stationary with re- 
spect to unitarian mixing among the orbitals. This lo- 
calization criterion is necessary, because the SIC is not 
invariant under unitary transformations of the occupied 
orbitals. This is in contrast with the LSD where a uni- 
tary transformation of the occupied orbitals leaves the 
LSD potential invariant, since the total charge density 
remains unaltered. For the orbital dependent SIC poten- 
tial V^^^ such a unitary transformation will change V^^^. 
The localization criterion < a|Vj^^^ - V^^'^IP >= de- 
termines the unitary transformation which ensures the 
global minimum of the total energy and the hermiticity 
of the Hamiltonian. Solutions of this equation usually 
take the form of the eigenvector \a > having weight in 
one channel only {{Imaj)'^ + (Reaj)'^ = 1) which would be 
different from the channel where the weight of the eigen- 
vector 1/3 > ((Im/3i)^ + (Re/3i)^ — 1) is concentrated, i.e., 
i is not equal to j. This generalization forms the ba- 
sis of the SIC implementationsi2*2i which start from a 
band-picture scenario. 



IV. SINGLE-SITE SIC-LSD FORMALISM 

As already mentioned, the proposed generalization of the 
Perdew and Zunger idea is based on the notion of reso- 
nances in scattering theory, which are the reminiscence 
of atomic states in the solid. Core states are represented 
by bound states at negative energies, where the imagi- 
nary part of the generalized complex phase shift jumps 
abruptly by n. Localized valence states still have very 
sharp resonances but band-like states are characterized 
by slowly varying phase shifts. 

The central quantity of (scalar-relativistic) multiple scat- 
tering theory is the single-particle Green's function— 

G.(r,r';e) = V Zi,(r,; e) r^^,(6) Zi,,(r' ; e) 



LL 



L 



Zl,(r<;e)JL(r>;6)<5. 



tj, 



(7) 



with r = Hi + Vi, where r.^ is a vector inside the cell at 
Hi, L = {l,m) denotes the combined index for the de- 
composition into symmetrized lattice harmonics Y^ and 
r<(r>) is the vector smaller (larger) in magnitude from 
the pair (r, r'). The building blocks of the Green's func- 
tion are the regular and irregular solutions of the radial 
Schrodinger equation at a given (complex) energy e 



ZlAr^■,e) = 


= ZUn;e)YL{ri) 


(8) 


zUn;^) = 


= ZUn-e)Yl{h) 


(9) 


JLin;^) -- 


= Jl{n-e)YL{h). 


(10) 



The scattering-path matrix t { in L,L' and a represen- 
tation) 



lie)=[r\e)-g{e)]- 



(11) 



is related to the structural Green's function g{e), describ- 
ing the free propagation between the scattering centers, 
and the t matrix defines the single site scattering. 
The total valence charge density per spin a is given by: 



can be written as 



1 /"^^ 

n„{r) = / deIniGo^(r,r;e), 

'^ Je„ 



(12) 



where Eb and Ep denote the bottom of the valence band 
and the Fermi energy, respectively. In standard LSD cal- 
culations, the new effective potential for the next itera- 
tion of the self-consistency cycle is calculated from this 
density (now including the core contributions) as 

Kff'?(r) = Kxt(r) + VH[n]{r) + V^J^^[nT,n^]{r). (13) 

In order to remove the spurious self-interaction, still 
present in this potential, we consider the problem of elec- 
trons moving in an array of scatterers. As already men- 
tioned, an electron which shows localized behavior has 
a sharp resonance in its phase shift, associated with a 
large Wigner-delay time on a particular site. To deter- 
mine the SIC charge we will consider for a moment the 
atomic limit, i.e., the situation where the scatterers are 
far apart. In this case the single site i-matrix and the 
local multiple scattering r-matrix coincide, and all occu- 
pied states correspond to bound states. In this limit each 
bound state contributes exactly the charge of one elec- 
tron, and this charge can be calculated by integrating the 
diagonal of the spectral function just around the energy 
of the bound state. In order to be able to decompose 
the charge density (Eq. 1)12(1 ') into different angular mo- 
mentum channels, we choose symmetry adapted spheri- 
cal harmonics. These are defined by applying a unitary 
transformation to the ordinary real (or complex) spher- 
ical harmonics, such that the on-site scattering matrix 
becomes diagonal. 



E 



uUrLLA^)UL,L'=SLL'fl\{e). (14) 



It is easy to verify that the required transformation ma- 
trix U is, in fact, independent of the energy e. This 
transformation to symmetry adapted spherical harmon- 
ics also ensures that the degeneracy of states, which are 
localized, is conserved. We will demonstrate this later 
by Sl-correcting the triplet states (one by one) of the 
Ce /-manifold. In this symmetrized representation, the 
Green's function, which in the atomic limit equals the 
atomic Green's function, becomes diagonal with respect 
to this quantum number. Hence we can decompose the 
spin resolved charge density into its L components and 
define the charge of a state, characterized by its principle 
quantum number n, angular momentum L and spin a: 



,sic 



1 



L;(r) = — 



deImGLL^cr(r, r;e) 



(15) 



where the energies Ei and E2 lie slightly below and above 
the energy of the state nLa. Within the multiple scatter- 
ing formulation, in the atomic limit, this charge density 



-^^(r) 



1 



-E2 



delm 






Zl,(r,;e)T^Je)ZL(r,;e) 



(16) 



where i is the site index, since in this case the single site 
t-matrix and the r-matrix are obviously identical. This, 
of course, is not the case for a solid with finite lattice 
spacings. When considering resonances in a solid it is 
not a priori clear whether to use the t- or r-matrix for 
calculating the charge density in question. The main dif- 
ference between using the t- matrix or r-matrix is that 
the latter does include a small hybridization of the local- 
ized state with the surrounding atoms, while the former 
does not. Also the choice of the lower and upper inte- 
gration limits is not clearly defined. We will now give a 
short discussion of the possible modes for calculating the 
SIC charge of a resonant state. 

The lower integration limit is most reasonably chosen to 
be the bottom of the energy contour E^. However, care 
has to be taken that this contour always encloses the SI- 
corrected states. Nota bene, in the case of Ce, discussed 
in the following sections, the contour also includes the 5p 
semi-core states. The upper integration limit could be 
either chosen such that the SIC charge density integrates 
to exactly one electron, or simply set to the Fermi energy 
E-p- Using the t- matrix, we find that we have to integrate 
up to an E'top, which is slightly above the Fermi level in 
order to capture one electron. The r-matrix, on the other 
hand, due to hybridization, yields a charge of one for 
energies slightly below the Fermi energy. Unfortunately, 
when dealing with the r-matrix, it is computationally 
very expensive to assume E2 different from the Fermi 
energy. However, the excess charge due to integrating up 
to the Fermi level is only of the order of a hundredth of an 
electron. (The missing charge in case of integrating the 
i- matrix up to the Fermi level is of similar magnitude.) 
In the following we used the r-matrix integrated up to 
Ey to determine the SIC charge. Some tests with the 
t-matrix, and the requirement of a SIC-charge of unity, 
resulted in an upward shift of the total energies by about 
1 mRy. 

The charge density, calculated in either of the proposed 
ways, is used to construct the effective self-interaction 
free potential, namely 



K 



SIC-LSD/ 
cGALa \^ 



(17) 
In this paper we only consider the spherically symmet- 
ric part of the SIC density and SIC potential. Hence, 
the i-matrix is diagonal in I and m. Here it should be 
noted that, if we transform the equations back to the un- 
symmetrized (real or complex) spherical harmonics, this 
effective potential assumes matrix character with respect 
to the angular momentum, and would not simply couple 
to the density, but rather to the non-diagonal Im, I'm' 



density matrix. This is conceptually analogous to the ro- 
tationally invariant formulation of LDA+U by Dudarev 
et alM 

For each self- interaction corrected channel L = (l, m) and 
CT, we replace the L — th element of the original t-matrix 
by the one obtained from the Sl-corrected potential 



tr,^ = tr,^(l - 5^l5a.d) + t'f^'^ ^^^5j^j6a,a, 



"-Lct 



^L,L^ 



(18) 



where t^^^. is the t-matrix calculated from the effective po- 
tential V^i^{v), and t^sic-LSD .^ calculated from the SI- 
corrected potential V^^^j^^ (r). This t-matrix is then 
used in Eq. l(TT|l to calculate the new, Sl-corrected, scat- 
tering path matrix r. From the latter the new SIC-LSD 
charge density is calculated, and the process is iterated 
until self-consistency is reached. The correction term, 
which approximately compensates the self-repulsion, is 
an attractive potential which will pull down in energy 
the state to which it is applied (see section rVTIll . 
To finish this section we would like to mention that in 
contrast to the LSD, the SIC-LSD Hamiltonian is not 
invariant under unitary transformations of the occupied 
orbitals. As pointed out before, in the full implementa- 
tion the localization criterion is applied to make the so- 
lution stationary under this unitary mixing of states. In 
the present implementation there is no such localization 
criterion, and one has to be solely guided by the ener- 
getics to find the global energy minimum. Note that the 
total energies are invariant under a rotation of the coor- 
dinate system owing to the symmetry adapted spherical 
harmonics that diagonalize the r-matrix at the F point. 
Hence, the energies of the configurations where, one by 
one, each state out of a degenerate manifold is localized 
are the same. This was tested on Ce by Sl-correcting all 
/-states separately. As expected, the energies for local- 
izing any of the Ti„ (or respectively the T2u) states were 
identical. 



V. CPA GENERALIZATION 

One of the advantages of the multiple scattering im- 
plementation of the SIC-LSD formalism is that it can 
be easily generalized to include the coherent potential 
approximation^2ii^4iiii^4i& extending the range of appli- 
cations to random alloys. In addition, one can use it to 
treat static correlations beyond LSD by studying pseu- 
doalloys whose constituents are composed e.g. of two 
different states of a given system: one delocalized, de- 
scribed by the LSD potential, and another localized, cor- 
responding to the SIC-LSD potential. Combined with 
the DLM formalism for spin-fluctuations, ■^■^'^^ this allows 
also for different orientations of the local moments of the 
constituents involved. 

In the CPA extension of the SIC-LSD formalism, bearing 
in mind its single-site aspect, it is required to satisfy the 



following CPA self-consistency condition 

cr^.oo(e) + (1 _ c)r^^o°(e) = r^^™(e), (19) 

where the impurity r— matrices T'^'°°(e) and T^''^°(e) are 
given by 



T^'™(e) = 



r^^™(e) 



(l + rC^0"(6)(i^(6)-tc(6))) 

Tg'0"(£) 

(l + r^'0"(e)(ts(e)-icW)) 



(20) 
, (21) 



and the r— matrix of the coherent potential approxima- 
tion 



T^' "(e) 



VL 



BZ 



d-^k 



{tc\e)-g{k,e)) 



(22) 



Here VI bz is the volume of the Brillouin zone (BZ), t^(e) 
and tg (e) are the respective single site scattering matri- 
ces of the A and B species, occurring with the concentra- 
tions c and 1 — c, respectively, and tp(e) is the i— matrix 
of the effective CPA medium. Note that in the CPA ex- 
tension of the SIC-LSD formalism, the CPA condition 
(Eq. H19|ll is an additional self-consistency criterion to 
the usual charge or potential self-consistency. 
Finally, it should be mentioned that the formalism of this 
section can be easily generalized from a binary to a multi- 
component case, as described in Ref. |23. In addition, an 
extention of the LSIC-CPA formalism to finite tempera- 
tures can be implemented as described in the next section 
(Ell). 



VI. FINITE TEMPERATURES 

In this section we summarize the relevant formulas under- 
lying the finite temperature generalization of the present 
formalism in its CPA extention. In contrast to T = 0, 
at finite temperatures the physics is dominated by ther- 
mal (classical) fluctuations. Therefore, to properly take 
into account the finite temperature effects, one needs to 
evaluate the free energy of the system (alloy) under con- 
sideration, namely 

F{T,c,V) = -Btot(T,c,y)-T(5ei(r,c,l/) 

+ 5'mix(c) -I- S'mag(c) -|- 5'vib(c) j , (23) 

where 5*01 is the electronic (particle-hole) entropy, S'mix 
the mixing entropy, 5'mag the magnetic entropy, and S^ih 
the entropy originating from the lattice vibrations. 
The electron-hole entropy is defined as^ 



5ei(T,c,y) = -ks Jden{e)(^ff,{e)\nfp{e) 

+ il-fpie))Hl-fp{e))), (24) 



where fcs is the usual Boltzmann constant and /^(e) de- 
notes the Fermi-Dirac distribution function. The entropy 
of mixing in the case of a binary system can be expressed 
as 



Sn.Uc) = -fcs(clnc + (1 - c) ln(l - c)). 



(25) 



The magnetic and vibrational entropies are strongly de- 
pendent on the system under consideration, and they will 
be discussed in more detail in the section dedicated to the 
phase diagram of Ce. 

Finally, note that in the definition of the free energy, 
the finite temperature enters only via the Fermi-Dirac 
distribution and the entropy contributions, while for the 
exchange-correlation energy, being part of the total en- 
ergy Etot7 the T=OK LDA (LSD) approximation is used 
for all temperatures, which is a common practice in all 
ab initio calculations. 

This section completes the formal description and im- 
plementation of the LSIC-KKR-CPA band structure 
method. In the following sections we shall illustrate the 
potential and versatility of this approach for describing 
strongly correlated electron systems by an application to 
Ce. There we present both the T=OK and finite temper- 
ature results, including the phase diagram of the famous 
a — 7 phase transition. 



VII. CE a - 7 PHASE TRANSITION 

Ce is the first element in the periodic table that contains 
an / electron, and shows an interesting phase diagrami^ii 
In particular, its isostructural (fee -^ fee) a — 7 phase 
transition is associated with a 15%- 17% volume col- 
lapse and quenching of the magnetic moment.'^'' The low- 
pressure 7-phase shows a local magnetic moment, and is 
associated with a trivalent configuration of Ce. At the 
temperatures in which the 7-phase is accessible, it is in 
a paramagnetic disordered local moment state. Increas- 
ing the pressure, the material first transforms into the a- 
phase, which is indicated to be in an intermediate valence 
state with quenched magnetic moment. At high pressures 
(50 kbar at room temperature) Ce eventually transforms 
into the tetravalent a'-phase. With increasing tempera- 
ture, the a — J phase transition shifts to higher pressures, 
ending in a critical point (600K, 20 kbar), above which 
there is a continuous crossover between the two phases. 
In the following paragraphs we first discuss the SIC and 
non-SIC /— phase shifts and densities of states. Then we 
compare the results of LSIC total energies, for the ferro- 
magnetic arrangement of local moments, with the earlier 
calculations of the full SIC implementation, in order to 
benchmark the new method. After discussing the density 
of states of the LDA and SIC-LSD calculations, we mix 
the two phases using the CPA and DLM for the spins. 
Finally, allowing for finite temperatures, we describe the 
full phase diagram of Ce. 
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FIG. 1; (color online) Phase shifts of the Sl-corrected and 
uncorrected /-states in Ce from the SIC-LSD calculation. The 
energies are relative to the Fermi level. 
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FIG. 2: (color online) Density of states and integrated density 
of states of the Sl-corrected /-channel. 



f-Phase shifts and corresponding densities of 
states 



Before presenting our results for the phase diagram, we 
discuss briefly the scattering properties of a single Ce 
site. In particular, we concentrate on the phase shifts 
and corresponding densities of states for / electrons. 
In Fig. ^we show the phase shifts for the Sl-corrected 
and uncorrected /-channels of Ce. It can be seen that 
the uncorrected /-states have a very sharp resonance just 
above the Fermi energy. The steep resonance corresponds 
to a long Wigner delay-time and indicates that the state 
is already well localized. The self-interaction corrected 
/-state is shifted down in energy by about 9 eV, and 
becomes a bound state (it lies below the muffin-tin zero). 
In Fig. 121 we present the density of states and the inte- 
grated DOS for the self-interaction corrected /-channel. 
It can be seen that by integrating Eq. (|16|) up to the Fermi 
energy, one collects slightly more than one electron. This 



is because there is a small density of states in the vicin- 
ity of the Fermi level (note the scale of the left-hand side 
axis in Fig. [21 , which is due to slight hybridization of the 
SIC channel with the other /-channels whose resonances 
occur in the vicinity of the Fermi level (see Fig.^. Some 
contribution to this density of states might also come 
from 5/ states. This is implied by the behavior of the 
phase shifts in Fig. ^ The sharp jump by tt indicates 
that 4/ state is a bound state of the SIC potential, and 
the slow rise of the phase shift above the Fermi energy 
can be associated with the progression towards the 5/ 
state. Figure |21 also shows that the integrated DOS at 
the energy where the phase shift goes through it, i.e., at 
about -2.5 eV, is slightly less than 1. This is most likely 
due to the integration method used to display the quan- 
tities in Fig. [21 which is less accurate than the contour 
integral used in the self-consistent calculations. 
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FIG. 3: (color online) The calculated total energies for Ce 
from LDA and SIC-LSD, with different /-states localized, as 
functions of volume, given in atomic units (a.u.)"^. 



In order to determine the ground state configuration of 
Ce at a given volume, we calculated the total energies 
for different volumes using the LDA to describe the a- 
phase and the SIC-LSD formalism for the 7-phase, when 
Sl-correcting one localized /-electron, allowed to popu- 
late in sequence all possible /-states. In both LDA and 
SIC-LSD calculations spin-orbit coupling has been ne- 
glected for valence electrons, but fully included for core 
electrons, for which the Dirac equation has been solved. 
The corresponding total energies as functions of volume 
are shown in Fig. [3 We find that the LDA, used to rep- 
resent the Qf-phase, yields the lowest energy minimum, 
as seen in Table Q] There the ground state properties 
of the all studied configurations are summarized. Ta- 
ble mi compares the present results for the ground state 
configurations with previous calculations and with exper- 
imental values. Note small differences between the differ- 
ent calculations, which are due to different schemes, and 
indicate the sensitivity of the results to computational 
details. 

The observed degeneracy of the states within the triplets 
demonstrates the rotational invariance of the formalism. 
Note the large crystal field (CF) splitting, separating the 
Tiu triplet from the other SIC-states. As already men- 
tioned, the calculations presented in this section assume 
a ferromagnetic alignment of the local moments in the 
7-phase. However, when discussing the phase diagram 
of Ce, we will also consider the disorder of the local mo- 
ments using the DLM frameworki22i24 
For the 7— phase, treated ferromagnetically, out of the 
three possible localized states listed in Table [3 the state 
with the A2U symmetry gives the lowest energy solution. 
This localized state is also associated with the highest 
volume among the all possible localized configurations. 
Only 0.8 niRy separate the minima of the a and 7 phases, 
giving rise to the transition pressure at the absolute zero 



TABLE I: The total energy differences as obtained from 
the LDA and SIC-LSD calculations, with respect to the 
ground state energy solution (LDA), for Ce in different /- 
configurations. The corresponding volumes and bulk moduli 
(evaluated at the theoretical lattice constants) are also given. 
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LDA 
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158 




701 


A2„ 


0.8 




202 




355 


SIC Ti„ 


20.3 




201 




352 


T2. 


1.5 




197 




351 



of about -2.3 kbar. This is in good agreement with the 
experimental value of -7 kbar, when extrapolated to zero 
temperature, and with other theoretical values (see also 
Table IIV|) . The bulk moduli, given in Table Q] are cal- 
culated at the theoretical equilibrium volumes. When 
evaluated at the experimental volumes (as it is common 
practice in DFT calculations), their values are substan- 
tially reduced to 239 kbar for the a-phase and 203 kbar 
for the 7-phase, which is in considerably better agreement 
with the experimental numbers. The volume collapse 
(with respect to the volume of the 7-phase) is obtained 
at 22%, which also compares well with the experimental 
values of 15%- 17%. We note that the underestimation of 
the volumes of both the a- and 7-phases is due to the 
KKR /—convergence problem, which was addressed by 
Moghadam et alm^ They demonstrated that angular mo- 
menta as high as 16 were needed to obtain satisfactory 
convergence in the total energy. In the present calcula- 
tions we chose /max — 3, which does not seem sufficient for 
a good description of the equilibrium volumes of the two 
phases. Although it seems that this /—convergence prob- 
lem should affect the LSD and SIC-LSD calculations in a 
similar manner, we see a significantly larger error for the 



TABLE II: Comparison of the computed equilibrium volumes 
and bulk moduli with those of other calculations and exper- 
iment. The bulk moduli have been calculated at the theo- 
retical equilibrium volumes. Note that unlike in the present, 
SIC-LSD (KKR) implementation, the results based on the 
LMTO refer to the full SIC-LSD scheme, involving repeated 
transformations between real and reciprocal spaces, and Bloch 
and Wannier representations. The two different sets of LMTO 
calculations refer to different basis sets and to different ways 
of solving the SIC-LSD eigenvalue problem. In the reported 
LDA and GGA calculations the 7-phase was modelled by con- 
straining the /-electrons to the core. 
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a-phase, in agreement with the results obtained by other 
well known KKR packages when the LDA approximation 
is implemented to describe the electronic structure of the 
a— phasei^ The larger error for the a-phase than for the 
7— phase (found also in the LMTO-ASA calculations) is 
most likely due to the fact that LDA is not adequate for 
describing the experimentally reported correlated nature 
of the a-phase. In fact, the LDA calculations correspond 
strictly to the high-pressure a'-phase, which is purely 
tetravalent and has a smaller lattice constant than the 
observed a-phase. However, as already mentioned, in 
our calculations we have treated Ce as a trivalent system 
(one localized /-electron) in the 7-phase, and a tetrava- 
lent system (all /-electrons are treated as delocalized) in 
the a-phase. Experimental data seems to suggest, that 
a-Ce has a non-integer valence of 3.67. One could argue 
that this intermediate valence character of the 4/ state 
could be represented in terms of a pseudoalloy composed 
of the trivalent and tetravalent Ce atoms. We shall elab- 
orate on this point in one of the following subsections. 



C. Densities of states of Ce a and 7 phases 

The densities of states of Ce from the LDA and ferro- 
magnetic SIC-LSD calculations are shown in Fig. 0] The 
LDA DOS shows all the /-states hybridized into the s-, 
p- and d-states. However, in the SIC-LSD panel of the 
figure, one clearly sees the split-off localized /-state at 
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FIG. 4: (color online) Spin-resolved density of states of Ce in 
the a- (a) and 7-phase (b) with ferromagnetic arrangement 
of local moments. 



about -8 eV. Of course, this docs not agree with the spec- 
troscopic position of this state. To accomplish the latter, 
one would have to take into account the self-energy, which 
could be evaluated from the total energy difference cor- 
responding to systems with constrained /-occupationsi^ 
The localized /-state apart, one can clearly see the ex- 
change splitting of the remaining states in the SIC-LSD 
calculations. The unoccupied /-states in the SIC-LSD 
density of states are pushed up by 1 eV or so. These 
unoccupied /-states are furthermore exchange split by 1 
eV. 

From Table IIIII we note that the number of occupied 
/-states hardly differs between the 7-phase and the a- 
phase. This non-integer value of 1.35 in both cases is a 
consequence of the hybridization of the /-states with the 
s, p and d states seen in Figs. and El The number 
of / electrons remains constant between the LDA and 
the SIC-LSD, which might seem rather surprising. What 
happens is that / electrons that participate in bonding 
in the a-phase get transferred to localized electrons in 
the 7-phase. At the same time some of the bonding 
d-electrons are transferred to the repulsive sp-channel. 
These effects conspire to give the larger lattice constant 
for the 7-phase. 

The s, p and d spin-resolved densities of states of Fig. [51 
show remarkably rigid band behavior between the LDA 
and SIC with the Fermi energy moving down with respect 
to the SIC-LSD d partial density of states. We also note 
that in Fig. the unoccupied /-states are well separated 
from the Fermi level, because the Fermi energy is lowered, 
and the hybridization with the occupied s, p and d states 
has been substantially reduced. These changes in the s, 
p and d densities of states are also reflected in Table IIIII 
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FIG. 5: (color online) Spin- and symmetry-resolved DOS 
(in states/eV) for the states originating from s{Aig)-, p{Tiu)- 
and d{T2g and £9)- channels. As in Fig. 21 the SIC-LSD cal- 
culation refers to the ferromagnetic arrangement of the local 
moments. 
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FIG. 6; (color online) Spin- and symmetry-resolved DOS (in 
states/eV) for the states originating from /- channels. Note 
the different scale of the plot with respect to Fig. |3 As in 
Fig. 0] the SIC-LSD calculation refers to the ferromagnetic 
arrangement of the local moments. 



TABLE III: Angular momentum decomposed charges from 
LDA and SIC-LSD calculations. Note that the p-channel in- 
cludes the 5p semi-core states. 



LDA 
SIC 



P 



f 



0.40 
0.51 



6.06 
6.16 



2.19 
1.99 



1.35 
1.35 



where we see a reduction of 0.2 electrons in the d-channel 
of the localized phase and a corresponding increase of 0.1 
electrons in both the s and p channels of the localized 
phase. Even though the Fermi energy moves down in the 
localized phase, we see from Fig. \^ that the number of 
occupied states in the s and p channels has increased. 



D. Ce a — 7 pseudoalloy 

In order to improve on the LDA representation of corre- 
lations in the a-phase of Ce, in the present approach, in 
the spirit of the Hubbard III approximation, 42, one can 
model the experimentally implied non-integer valence of 
the Ce ions by a pseudoalloy consisting of the trivalent 
(SIC-LSD) Ce ions with concentration c, and the tetrava- 
lent (LDA) Ce ions with the concentration (1-c). In ad- 
dition, taking into account the disordered local moments 
of the trivalent Ce ions in the 7-phase, one can assume 
that their up and down orientations occur with equal 
probabilities. Supposing homogeneous randomness, such 
a ternary pseudoalloy can be described by the coherent 
potential approximation (CPA). The respective concen- 
trations of the trivalent and tetravalent Ce ions in the 
pseudoalloy are then determined by minimizing the total 
energy for each volume with respect to the concentration 
c. 

In Fig. [3 wc show the total energies for the a-j pseu- 
doalloy at T=GK, in which the 7-phase occurs with the 
concentration c (c/2 for each spin orientation), and the 
a-phase with the concentration (1-c), for several lattice 
constants. It can be seen that all shown total energy 
curves have their minima either at c=0 (pure a-phase) 
or c=l (pure 7-phase). Hence a fractional occupation 
of the 4/ state appears to be energetically unfavorable 
for all lattice constants. From these calculations we can 
conclude that a static, single-site approximation is not 
sufficient to describe the intermediate valence state of a- 
Ce at T=0 K. These calculations are consistent with the 
earlier results by Svane who performed supercell calcula- 
tions to model 25%, 50% and 75% of a-7 admixtures, but 
treating the 7-phase ferromagnetically and not as a DLM 
phasc,"^^ There too, no total energy minimum was found 
for intermediate concentrations between and 1, and also 
a mainly convex (from above) curvature for the total en- 
ergy, as a function of concentration, was obtained. This 
suggests that to describe the intermediate valence state of 
the a-phase one would need to consider a dynamical gen- 
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eralization of the CPA^?^ which would involve dynamical 
fluctuations between the trivalent and tetravalent states. 
Other possible mechanisms to favor intermediate valence 
will be commented about in section rvilll 



E. Ce a — 7 spectral functions 

In this section we discuss the spectral functions along the 
V — X line, calculated for the pseudoalloy consisting of 
50% 7-admixture into the a-phase, i.e., 50% a, 25% 7 
spin-up and 25 % 7 spin-down, in comparison with the 
spectral functions of the pure phases, all at the same vol- 
ume, as shown in Fig. |S1 The pure 7-phase has been 
represented by a 50 % spin-up and 50 % spin-down al- 
loy. The pure a-phase (panel (a)) shows a well defined 
band-structure. (The minor smearing of the bands is 
due to a small imaginary part added to the energy). Of 
course, the LDA leads trivially to a non spin-polarized 
band structure, but the absence of an exchange splitting 
in the panels (b) and (c) is due to the use of DLM, which 
defines an effective medium in which the local moments 
are averaged out. In the panels (b) and (c) the broad- 
ening of the spectral functions is apparent. The actual 
line- width of the spectral function can clearly be seen in 
Fig. |5| showing the spectral functions at the F point. 
Here, similarly to the smearing effect seen in Fig. |H1 the 
residual linewidth seen in panel (a) is purely due to the 
small imaginary part of the energy, necessary to obtain 
well behaved r-matrices. 

The spectral functions reveal some features of the cur- 
rent approach. Firstly, it can be seen that the s- and 
d-derived states at the F-point are hardly affected by 
the CPA. These states have no hybridization with the 
Sl-corrected /-state. Symmetry analysis of the spectral 
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function, shown in panel (b) of Fig. ^ reveals that the 
^2m /-state appears twice (the peak at -8 eV and the 
sharp shoulder just above the Fermi energy), since we 
are in the split-band regime. Due to the DLM treat- 
ment of the 7-phase, this feature is also seen in panel (c), 
where the upper of the split-band peaks merges with the 
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FIG. 9: (color online) Spectral functions at the F point for the 
7-phase concentrations (a), 0.5 (b) and 1 (c). The character 
of the respective peaks is marked above the upper horizontal 
line of each panel. 



lower triplet /-peak. The two triplets {Tiu and T2„) show 
common band behavior. The corresponding broadenings 
are very different between panels (b) and (c). Moreover, 
panel (c), as compared to panel (a), also shows that the 
unoccupied /-states have been pushed up in energy, be- 
cause the localized /-electron is more effective in screen- 



ing the nuclear charge. This results in an energy split- 
ting of the a and 7 unoccupied triplets which leads to a 
broadening of the unoccupied triplets, as seen in panel 
(b). The shoulder at about 1 eV is an indication that 
the splitting is noticeable on the scale of the dispersion 
of the bands. The broadening of the triplets is reduced 
in panel (c) where no a-phase is admixed and where the 
broadening of the lower triplet state is a consequence of 
merging with the upper split-band peak of the singlet 
state. Note that the /-states at the Fermi-level have a 
finite lifetime which might indicate a shortcoming of the 
static CPA for the description of an intermediate valence: 
a coherent mixture of the localized and delocalized states 
in terms of wave functions might be a more appropriate 
description. 



F. Phase diagram 

In this subsection we concentrate on the finite temper- 
ature phase diagram of Ce. The idea of describing Ce 
at finite temperatures as a pseudoalloy of a- and 7-Ce 
atoms was first put forward by Johansson et alm^ and by 
SvaneH^ Since at finite temperatures the thermal (classi- 
cal) fluctuations are of major importance, the static ap- 
proximation should suffice. In the work by Johansson et 
al. the pseudoalloy was treated by the CPA implemented 
within the LMTO method, where the 7-phase was mod- 
eled by including the 4/-states into the core, while in the 
a-phase the / states were treated as band-states. Due 
to the different treatment of both phases, their total en- 
ergies could not be compared and the energies of the 
7-phase had to be adjusted by hand to yield the correct 
a — 7 transition pressure at zero temperature. Svane, on 
the other hand, described the 7-phase as a ferromagnet 
using the LMTO-SIC, thus treating both phases on equal 
footing, and utilizing a supercell to mimic the pseudoal- 
loy at only a few accessible concentrations. From these 
calculations he concluded that a linear interpolation of 
the a- and 7-energies to arbitrary concentrations should 
be adequate enough. 

Here we present calculations which combine both the 
CPA and the SIC-LSD to describe the 7-phase as a DLM 
system, treated as a ternary alloy, consisting of spin-up 
and spin-down SIC sites with concentrations c/2 each, 
and LDA sites with the concentration (1 — c). In addi- 
tion, we go beyond the scope of previous works by taking 
into account the effect of finite temperatures on the elec- 
tronic total energies and the electronic contribution to 
the entropy, as defined in section lvTI However, the vibra- 
tional entropy, Svib, is neglected in the presented results. 
We shall briefly comment on its effect on the calculated 
phase diagram in the next subsection where we analyze 
in detail how different aspects of the present calculations 
influence the final results. 

Ideally, one would like to treat a pseudoalloy, which con- 
sists of all possible states of a Ce-ion, i.e. the LDA, and 
all possible Si-corrected / states. This would give rise 
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to a pseudoalloy consisting of 15 components: the LDA 
plus the 14 possible /-SIC states (including the spin mul- 
tiplicity). Since this would be a quite formidable task, we 
use a simplified approach. FigureOlindicates that crystal- 
field splitting gives rise to nearly degenerate A2u (singlet) 
and T2u (triplet) states, while the Tiu triplet lies 20 mRy 
higher in energy. At the temperatures considered here, 
the Tiu states are thermally not accessible. Thus treating 
the remaining eight states as degenerate, leaves us with 
a nine component pseudoalloy, with the constraint that 
the concentrations of the considered eight SIC states are 
equal and can be set to c/8. It is easy to show that in 
this case, in addition to the mixing entropy defined in 
section IVll one has to take into account a term 



r=OK 



5(c) = fcsclnS, 



(26) 



arising from these eight-fold multiplicity. Note that in 
previous studies, where the CF splitting has not been 
taken into account, the magnetic entropy was assumed 
to be that of a spin-orbit (SO) coupled J = 5/2 state, i.e. 
»S'mag(c) — /csclnG. In the next subsection we shall com- 
ment on how the two different magnetic entropy terms in- 
fluence the critical characteristics of the calculated phase 
diagram. 

We performed calculations for several lattice constants, 
embracing the equilibrium lattice constants of both 
phases, and concentrations from to 1, in steps of 0.1. 
The results for the free energies of the three selected tem- 
peratures are shown in Fig. ^| In the T = OK panel one 
clearly sees the two minima, corresponding to the pure 
a— phase (LDA) and the pure 7— phase (SIC-LSD) calcu- 
lations. The equilibrium lattice constant for a given con- 
centration c interpolates between the two extremes, and 
it is apparent that an intermediate valence state, even if 
not energetically favorable, would correct the underesti- 
mated lattice constant of the a-phase. As the temper- 
ature is increased, the free energy surface gets strongly 
tilted towards the SIC-LSD side and now shows only one 
broad minimum. This is mainly due to the magnetic 
entropy. 

Although the concentration-dependent free energies are 
the quantities directly accessible from the calculations, 
they do not easily reveal the full information on the phase 
diagram. In order to determine the full p — T phase dia- 
gram, it is necessary to calculate the Gibbs free energy: 



GiT,c,p)^FiT,c,V{T,c,p))+pViT,c,p). (27) 

The Gibbs free energies are displayed in Fig. ^2 From 
them, at each given pressure and temperature, we can 
determine the concentration of the trivalent Ce, by min- 
imizing the Gibbs free energy with respect to c. At zero 
temperature one finds (for low pressures) two local min- 
ima, associated with c = and c = 1. By increasing the 
pressure, the order of the minima changes and the mini- 
mizing concentration jumps from 1 to 0. At higher tem- 
peratures the minima start moving towards intermediate 
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FIG. 10: (color) Calculated free energies for the temperatures 
r=0, 800 and 1600K. The plots represent fits to the calcula- 
tions, which have been performed for concentrations from 
(corresponding to the pure a-phase) to 1 (corresponding to 
the pure 7-phase), in steps of 0.1, and for lattice constants 
from 8.25 a.u. to 9.65 a.u., in steps of 0.2 a.u. A constant of 
17717 Ry has been added to all energies. 
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FIG. 11: (color) Gibbs free energies for T=0, 800 and 1600K. 
In order to enhance the readabihty of the plots the energies 
have been calibrated by a linear term, proportional to pres- 
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FIG. 12: (color online) The free energies as function of the 
volume for the temperatures (highest curve), 200, 400, 600, 
800, 1000, 1200, 1400 and 1600 K (lowest curve). A constant 
of 17717 Ry has been added to all energies. 



concentrations. Only above the critical temperature, one 
finds the niinimum smoothly changing from low to high 
concentrations . 

We can obtain the free energy of the physical system at a 
given volume by evaluating the concentration dependent 
free energy at the minimizing concentration Cmin: 



F{T, V) = F{T, 



.V). 



(28) 



These free energies are displayed in Fig.^J which clearly 
shows the double-well behavior for low temperatures, 
which is gradually smoothened out with increasing tem- 
peratures. Furthermore one finds that, at elevated tem- 
peratures, the free energy is mainly lowered at large lat- 
tice constants, corresponding to the 7-phase, with its 
larger entropy. 
Inserting the minimizing concentration c^in into the 
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FIG. 13: (color) Calculated isotherms for the temperatures 
T=0 (lowest curve), 200, 400, 600, 800, 1000, 1200, 1400, 
and, 1600K (highest curve). The color indicates the fraction 
of localized electrons: blue is all localized (7-phase) and red 
is all delocalized (a-phase). 
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pressure- volume relation 

p{T, V) = p{T, cmin, V) = -^F{T, c,„in, V) (29) 

allows to calculate the isotherms of Ce, which are dis- 
played in Fig. El It can be seen that the average valence, 
close to the coexistence line, gradually changes with in- 
creasing temperature. Above the critical temperature, 
the valence changes continuously with increasing pres- 
sure from trivalent to tetravalent. 
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FIG. 14: (color) Phase diagram obtained for the pseudoalloy, 
composed of a- and 7-Ce. The crosses indicate the calculated 
and experimental critical points. The color indicates the frac- 
tion of localized electrons: blue is all localized (7-phase) and 
red is all delocalized (a-phase). 

In Fig. 1141 we present the phase diagram, obtained from 
the free energies of the a-7 pseudoalloy, with the 7-phase 
described by the DLM approach. It can clearly be seen in 
the figure how the transition becomes continuous above 
the critical temperature. The experimentally observed 
critical point (600K, 20 kbar) falls on top of the cal- 
culated phase separation line, which starts at the zero 
temperature transition pressure of -7.4 kbar. This means 
that the slope of the phase separation line is in very good 
agreement with experiments. The critical temperature 
overestimates the experimental one by roughly a factor 
of 2, which is still reasonable considering that the critical 
temperature is very sensitive to various small details of 
the calculations and in particular the theoretical lattice 
parameters of both the Ce phases. Note that the Tc at 
zero pressure of 169 K (see Table HVll compares well with 
the experimental value of 141±10 K. 
Finally we examine in more detail the discontinuity 
across the phase separation line. Figure El shows the 
magnitude of the discontinuities for the various ingredi- 
ents of the Gibbs free energy. As expected, all contri- 
butions vanish at the critical temperature, above which 
there is a continuous cross-over between the a- and the 
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FIG. 15: Discontinuities of the total energy (thick solid line), 
the total entropy TS (thick dashed line) and the pV term 
(thick dashed-dotted line) over the phase separation line as 
function of the temperature. The entropy term is further 
decomposed into the electronic (thin solid line), the mixing 
(thin dashed line) and the magnetic (thin dashed-dotted line) 
contribution. 



7-phase. It also can be seen from this figure that the en- 
tropy discontinuity is by far the largest contribution. The 
phase transition is therefore driven by entropy, rather 
than by energetics. The entropy discontinuity itself is 
mainly determined by the magnetic entropy. 



G. Analysis of results 

In order to investigate the importance of the different as- 
pects of the present calculations (i.e. the DLM descrip- 
tion of the 7-phase, the inclusion of finite temperature 
effects in the electronic free energy, and the CPA itself), 
as compared to earlier studies, we have also performed a 
set of calculations where, selectively, we neglect some of 
these effects and look at the consequences. In particular, 
we study the influence of these effects on the critical tem- 
perature and the slope of the phase separation line. The 
results of these calculations are summarized in Table Hvl 
in comparison with the results of earlier theoretical, as 
well as, experimental studies. 

First we focus on the importance of the disordered lo- 
cal moments representation of the 7-phase. By compar- 
ing the DLM-results in Table IIVI with those marked by 
'Ferro' (ferromagnetic calculations for the 7-phase), one 
finds that the DLM calculations lead to a moderate low- 
ering of the critical temperature, and a more negative 
zero temperature transition pressure. This can easily be 
understood, since at zero temperature, the ferromagneti- 
cally ordered phase has a lower energy as compared to the 
disordered phase. Experimentally such a magnetic order 
is not observed, since at low temperatures (and positive 
pressures) Ce is in its non-magnetic a-phase. The low- 
ering of the critical temperature cannot be easily identi- 
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TABLE IV: The critical temperature and pressure, as well as, the zero temperature and room temperature transition pressures 
and the zero pressure transition temperature for different calculations: DLM and Ferro refer to the disordered or ferromagnetic 
alignment of the local moments in the 7-phase, CF and SO indicate the crystal field or spin-orbit scenario, as discussed in 
the text. The index I denotes calculations with finite temperature effects included in the band-structure and the CPA, II 
refers to the neglect of these finite temperature effects, and III represents calculations, where in addition the concentration 
dependence was approximated by a linear interpolation. The main results, which are also shown in the figures, are the DLM-CF 
(I) calculations, printed as bold in the table. 
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fled with a specific aspect of the DLM calculations, since 
many efi'ects, such as the curvature of the free energies 
with respect to the concentration, but also the anhar- 
monic terms in the total energy as a function of volume, 
conspire to determine the phase diagram. 

The second point is the effect of finite temperature on 
the total energy and the electronic entropy. In ear- 
lier studies of the phase diagram;22iSiiii the electronic 
structure calculations were performed at zero temper- 
ature, and the finite temperature effect entered only 
via the mixing entropy, the magnetic entropy and in 
Ref. I33 also the vibrational entropy. This means that 
Etot{T,c,V) was replaced by £'tot(0, c, V^) and the elec- 
tronic entropy Sd was neglected. In Fig. 1161 we an- 
alyze the difference between the electronic free energy 
Fei(r, c, V) = Etot{T, c, V) - TS'oi(r, c, V) and the total 
energy at T = 0. The difference Fei{T, c, V) ~ Fci(0, c, V) 
exhibits a moderate dependence on the concentration. 
The larger effect for the a-phase is easily explained us- 
ing a low temperature expansion. The main effect of 
the finite temperatures is the broadening of the Fermi 
function. To lowest order in temperature, the change of 
the free energy is proportional to the density of states 
at the Fermi level. The effect on the phase diagram can 
be seen in Table Hvl by comparing the columns I and II. 
Neglecting these finite temperature effects gives rise to 
an increase of the critical temperature by roughly 200K, 
while the slope of the phase separation line remains un- 
altered. 

Next is the effect of the CPA. In the study by Svane,'^^ 
it was suggested that the weak departure from linearity 
of the total energy curves as a function of concentra- 
tion (Fig. 0) did not play an important role. Thus, we 
recalculated the phase diagram, replacing the full con- 
centration dependence of the total energy by the linear 
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FIG. 16: Temperature dependence of the electronic free en- 
ergy as defined in the text. 



interpolation 



Sun(c, V) = {1- c)E{0, V) + cE{l, V). 



(30) 



The effect of this approximation can be seen in Table 
IIVI by comparing the columns II and III, respectively. 
Both sets of calculations use the T — OK total energies 
only. The transition temperature obtained from the lin- 
ear interpolation is strikingly reduced, in comparison to 
the CPA calculation, and similar to the one reported by 
Svane. As pointed out by Johansson et al.,^^ the criti- 
cal temperature is to a large extent determined by the 
mixing entropy. Without the mixing entropy one would 
at all temperatures find a sudden transition (with a fi- 
nite volume collapse) between the low and high pressure 
phases. The mixing entropy will, if the temperature is 
high enough, lead to a minimum of the free energy for 
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an intermediate concentration of the 7-phase, eventu- 
ally resulting in a continuous crossover between the low- 
pressure and the high pressure phase. As seen in Fig. 
the CPA gives rise to a convex (from above) curvature 
of the total energies. It effectively reduces the mixing 
entropy and increases the critical temperature. It can 
also be seen that within this linear approximation, the 
difference between the DLM and the 'Ferro' calculations 
is strongly reduced, which means that the main effect of 
DLM is not only the energy lowering of the 7-phase, but 
more importantly a different shape of the energy versus 
concentration curves. 

Last, we discuss the effect of the degeneracy of the 7- 
phase ground state. As described above, in the CF cal- 
culations a degeneracy of 8 was used, corresponding to 
almost degenerate A2U and T2U states, arising from a 
crystal field splitting. The results to be compared here 
are shown in the columns CF and SO of Table Hvl In the 
SO scenario, the term proportional to In 6, corresponding 
to a J = 5/2 state, was used for the magnetic entropy, in- 
stead of In 8, as in the CF scenario. One can see that the 
magnetic entropy, being linear in c, determines mostly 
the slope of the phase separation line and has only a mi- 
nor effect on the magnitude of the critical temperature. 
The results obtained with the DLM and the magnetic 
entropy due to the CF splitting show better agreement 
of the calculated transition pressures with experiment. 
This, however, may be due to a cancellation of errors. 
The slope of the phase separation line is mainly given by 
the ratio (S^ — Sa)/iVj — Va), where Sa:(j) and Vq(^) are 
the the T=0 values of the entropy and the volume of the 
a(7)-phase. 

Due to the larger underestimation of the volume of the 
a-phase, the volume collapse is slightly overestimated in 
our approach. Therefore the higher value of the entropy 
difference in the CF scenario, as opposed to the SO sce- 
nario, leads to a better agreement with the experimental 
slope. It should be noted though that these calculations 
do not include the vibrational entropy. To estimate the 
effect of this vibrational entropy, we recalculate the phase 
diagram using a simple model for the vibrational entropy, 
namely 



2000 



5'vib(c) = -fcscAS'Jii, 



'7 — a 



(31) 



with the value of AS'Jj^" « 0.75, as suggested by Jeong et 
al.^^ In doing so, the critical temperature is only slightly 
reduced to 1292 K, while the critical pressure increases to 
69 kbar. In view of the above discussion, it is not surpris- 
ing that a contribution, which is purely linear in c, basi- 
cally affects only the critical pressure. We conclude from 
this result that this simplified model for the vibrational 
entropy is too crude and a more sophisticated one should 
be put in place. One should here add that Johansson 
et alm^ have employed a Debye-Griineisen model, while 
Svane^ has completely neglected the vibrational entropy. 
Finally, we would like to comment on the finding of 
Niklasson et al.'^^ that the disordered local moments can 
give rise to a localization of f-electrons. From our DLM 
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FIG. 17: Critical points of the phase transition, obtained 
when the concentration of trivalent Ce atoms in the Q-phase is 
artificially fixed at a finite cq. The points are marked by their 
corresponding value of co, where co = represents purely 
tetravalent Ce. 



only calculation, without applying SIC, we found that 
it is not possible to stabilize a local moment in Ce, ex- 
cept at very large volumes, where even the LSD yields a 
magnetic solution. 



VIII. DISCUSSION 

As already mentioned, experiments indicate that the a- 
phase is not composed of tetravalent Ce atoms, but is 
rather described by an intermediate valence of 3.67i^ 
From Fig. 1101 we see that an intermediate valence, i.e., 
a fractional concentration (0 < c < 1), would lead to an 
increased equilibrium lattice constant. One can simulate 
the effect of the intermediate valence for the a-phase by 
simply rescaling the concentrations when evaluating the 
phase diagram. In Fig. El we present the critical points 
of the phase transition obtained when we represent the 
a-phase by a non-zero concentration, cq, of the trivalent 
Ce atoms, added into the host of tetravalent Ce atoms. 
As can be seen in the figure, the critical temperature 
quickly decreases with the increase of the admixture of 
trivalent Ce in the a-phase. In this calculation, the SIC- 
LSD energies have been uniformly calibrated to keep the 
zero temperature transition pressure at its original value 
of -7.4 kbar. It can be seen in the figure that the best 
value for critical temperature is obtained for cq ~ 0.4, 
which corresponds to an intermediate valence of 3.6. 
The intermediate valence scenario for the a-phase, as dis- 
cussed above, could result from dynamical fluctuations. 
These fluctuations could be realized by describing Ce as 
a two level system (TLS). We will elaborate on this idea 
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in section IIXI A mechanism, based on the dynamical 
interaction of two states will, quite generally, be more 
effective if the two states are close in energy. Looking 
at the total energy as a function of the lattice spacing 
and the concentration fFig. IIOJI . one finds that the pure 
a- and 7-solutions (LDA and SIC-LSD, respectively) are 
degenerate close to a = 8.9 atomic units. The interac- 
tion of these two states will lower the total energy at 
these lattice constants, and eventually might establish 
them as the global minimum. The increased lattice con- 
stant would be closer to experimental values. Clearly 
such a state would be neither described by full localiza- 
tion nor delocalization of the /-electron, but would be 
better characterized as an intermediate valence state. In 
the following section we will outline how we envision this 
two-level system to work. 

The next point we want to discuss is the effect of lat- 
tice relaxations. In the intermediate valence regime (as 
described by the static CPA) there is a rather large size- 
mismatch between the a- and the 7-atoms, which may 
give rise to strong internal strains. Allowing for lattice 
relaxations, which are not considered here, would give 
rise to an energy lowering for intermediate concentra- 
tions and could lead to intermediate valence, even in the 
static limit, and hence to the reduction of the critical 
temperature. 

Another factor that might have a significant influence on 
the phase diagram and its characteristics is associated 
with the single-site aspect of the CPA. Being a single- 
site theory, the CPA cannot deal with 'order in disorder', 
namely with short range order (SRO) in the distribution 
of the 'alloy' configurations. However, even in the dis- 
ordered phase it must be important to distinguish sit- 
uations where the nearest neighbours are preferentially 
like atoms from those where they are unlike atoms. Such 
short range order is expected to lower the free energy of 
the disordered state and hence lower Tc of the a— 7 phase 
transition. In addition, it should influence the energeti- 
cally favourable relative concentration of the two different 
components of the alloy. Thus by taking SRO into ac- 
count one might be able to move beyond the primitive al- 
loy analogy and improve on the present LSIC-KKR-CPA 
approach. The way it could be accomplished is by imple- 
menting the nonlocal extention of the KKR-CPAj^^i^Si^ 
which would allow to treat possible correlated valence 
fluctuations near the a ^ ^ phase transition in Ce. 

Finally, the phase diagram of Ce suggests that at nega- 
tive pressure there could exist a quantum critical point 
(QCP), i.e., a localization - delocalization transition at 
zero temperature, driven by pressure. The vicinity of this 
QCP, although it is not accessible experimentally, could 
still influence the physics of the material in the accessi- 
ble positive pressure range. The quantum fluctuations, 
which are responsible for the transition, should also be 
visible in its vicinity and could explain the correlated 
nature of the a- phase. 



IX. OUTLOOK 

Several times in this article we have referred to the pos- 
sibility of going beyond the theory outlined above by al- 
lowing for dynamical valence and spin fluctuations. In 
short, the suggestion is that, as above, we regard the 
self-interaction corrected and the not so corrected ver- 
sion of the local potential at a site, as corresponding to 
two states of the 'atom' and allow for such 'atom' to tun- 
nel between the two states \a > and |& >, which would 
form a two-level system. Electrons interacting with such 
TLS were already thoroughly studied in the context of 
metallic glassesi^ In the present scenario such a pro- 
cedure would address the valence fluctuations. On the 
other hand, identifying the possible spin states of the SI- 
corrected system with the two levels^ would constitute a 
possible dynamical generalization of the DLM formalism, 
taking into account dynamical spin fluctuations. This 
might describe the Kondo screening of the local moments 
at low temperatures. 

Many of the consequences of such interactions are by now 
well known44 In such studies the TLS is an atom tunnel- 
ing between two nearly degenerate positions in a metallic 
environment.— The physics is particularly interesting in 
the case of 'assisted tunneling' where the TLS changes its 
state as an electron scatters from it. However, it should 
be stressed that in the cases where the atom changes its 
position in the tunneling process, the TLS is external to 
the electron system whilst in the proposed model it is 
designed to capture the physics of a slowly changing col- 
lective degree of freedom of the electrons themselves. In 
this sense our TLS is very much like the DLM in para- 
magnetic metals. 

From the point of view of the present perspective the 
most directly relevant work on external TLS's is that 
of Vladar et alM- They show that, if the TLS couples 
two or more angular momentum states, as the tem- 
perature is lowered the high temperature 'pseudoalloy' 
phase gives way, at a characteristic temperature T^^-, to 
a ground state which is an orbital analogue of Kondo 
singlets of magnetic impurities in metals. In the case of 
the LDA/SIC-LSD TLS, this means, that at high tem- 
peratures, we would flnd the pseudoalloy, as described 
above, while at low temperatures, the system would be 
neither in the LDA, nor in the SIC-LSD state, but would 
be in an intermediate valence ground state. Clearly, such 
dynamical fluctuations may also help to reduce the crit- 
ical temperature if the generalized Kondo-temperature, 
associated with the TLS, is not too small. 
As a final remark concerning the dynamical generaliza- 
tion of our disordered local valence calculations we note 
that the above model is an analogue of Yuval and An- 
derson's 'Kondo Hamiltonian' approach to the magnetic 
impurity problem, as opposed to the full dynamical cal- 
culation based on the Anderson model, deployed for the 
same problem by Hamann42iffi The relevant point to 
stress is that these two calculations yield, in the ap- 
propriately asymptotic, namely scaling, regime the same 
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results. Thus, they lend credit to the above proposed 
short-cut to a first principles DMFT treatment of our 
fluctuating valence problem. 



X. CONCLUSIONS 

We have presented a new, multiple scattering, implemen- 
tation of the SIC-LSD formalism for solids within KKR 
band structure method, combined with the CPA descrip- 
tion of intermediate valences. The method has been il- 
lustrated on the application to the Ce a — 7 phase tran- 
sition. The results have been discussed in detail, high- 
lighting the new functionality and potential of this ap- 
proach owing to a better static description of spin and 
valence fluctuations. The importance of all the different 
aspects of the formalism has been analysed in detail. This 
new method is not to be looked at as an alternative to 
the earlier implementations within the LMTO-ASA band 
structure method. Its great potential, and in some way 
superiority, arises from the local and multiple scatter- 
ing aspects through which the method lends itself easily 
to various generalizations and extentions on the account 
of the straightforward determination of the one-electron 
Green's function. Of particular interest here is an inclu- 
sion of dynamical fluctuations, for which a roadmap has 
been briefly outlined in the previous section. The results 
of the present paper constitute the first crucial steps on 
this road towards dynamics. 
Finally, there is one more aspect of the present results 



which warrants further comment. As we have reported in 
Sec. lVIim in our DLM calculations local moment formed 
only when the local electronic structure was described 
by SIC-LSD, whilst it iterated to zero when LDA was 
used in recalculating the spin polarized crystal potential. 
This is precisely the behavior one would have expected on 
the basis of the numerous successes of SIC in predicting 
moment formation, and no moment formation in appli- 
cations to extended systems i^ii^S Clearly, the fact that 
our local implementation of SIC, namely LSIC, behaves 
in this way lends strong support to our contention that 
self-interaction correction formalism can be, and perhaps 
should be, applied at the local level. 
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